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This paper describes and compares two methods for estimat- 
ing the variance function associated with iTRAQ (isobaric tag for 
relative and absolute quantitation) isotopic labeling in quantitative 
mass spectrometry based proteomics. Measurements generated by the 
mass spectrometer are proportional to the concentration of peptides 
present in the biological sample. However, the iTRAQ reporter sig- 
yy ■ nals are subject to errors that depend on the peptide amounts. The 

C/3 ' variance function of the errors is therefore an essential parameter for 

evaluating the results, but estimating it is complicated, as the number 
of nuisance parameters increases with sample size while the number 
^ I of replicates for each peptide remains small. Two experiments that 

^f^ . were conducted with the sole goal of estimating the variance function 

(^ ' and its stability over time are analyzed, and the resulting estimated 

^^ , variance function is used to analyze an experiment targeting aber- 

^T ' rant signaling cascades in cells harboring distinct oncogenic muta- 

^^ ' tions. Methods for constructing conservative p- values and confidence 

C~^ , intervals are discussed. 

cn 

1. Introduction. Improvements in mass spectrometer resolution, accu- 
racy and sensitivity coupled with the development of increasingly sophisti- 
cated algorithms for protein identification from spectra have resulted in mass 
spectrometry (MS) becoming the tool of choice in large scale proteomics re- 
search. Typically, the mass spectrometer is used to measure short portions 
of the proteins called peptides. These are subjected to a process called tan- 
dem mass spectrometry (MS/MS) which ultimately yields a mass spectrum 
containing peaks which correspond to the primary amino acid sequence and 
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Fig. 1. Left: scatterplots of Yn versus Yfz m the March experiment. Right: pomtwise 
estimates of mean and variance and different models for the variance. 



enable the identification of peptides. When samples are labeled with iTRAQ 
(isobaric tag for relative and absolute quantitation) stable isotope reagents, 
the MS/MS spectra also contain peaks at predefined masses, whose intensi- 
ties provide a relative measure of the peptide abundance in a set of samples. 
A single experiment can yield tens of thousands of spectra identifying thou- 
sands of peptides belonging to thousands of proteins. Analysis of samples 
from different sources, for example, cells expressing different mutations in 
a known oncogene versus the wild- type (e.g., nontransforming) counterpart 
can provide insight as to the mechanisms by which specific genetic lesions 
manifested in the same protein drive a malignant phenotype. A workflow 
diagram and a brief explanation of the iTRAQ technique are given in Part 
A of the supplemental article [Mandel et al. (2013)]; for a detailed discussion 
on MS techniques see Eckel-Passow et al. (2009). 

In an experiment conducted in March 2009 and described below, Zhang 
et al. (2010) applied the iTRAQ protocol using two different labels for the 
same biological sample. The experiment yielded 2174 pairs of measurements 
corresponding to the amounts of 2174 peptides in the sample. The left panel 
of Figure 1 depicts the data^ (on a logarithmic scale) and clearly shows that 
the variance of peak intensity measurements is nonconstant and depends on 
the mean. This variance should be estimated for better understanding of 
the results of the MS analysis and to enable statistical inference about the 
peptide amounts. The right panel of Figure 1 displays the mean of each pair 
against its variance together with several estimates of the variance function 
described in Section 2.3. 



^Pairs with exactly the same reahzed value were excluded from the analysis; there were 
two such peptides. 
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The purpose of the current paper is twofold: first, to estimate the vari- 
ance function related to observations obtained by the common technique of 
iTRAQ stable isotope labeling [Ross et al. (2004), Aggarwal, Choe and Lee 
(2006)], which provides measurements of the relative amounts of peptides 
from two different biological samples in a single experimental run; and sec- 
ond, to construct confidence intervals for the abundance of a given peptide 
in a biological sample and to the ratio of two abundances under different 
conditions (e.g., cancer and wild type cells) using the estimated variance 
function, and to calculate p-values for the hypothesis of equality in peptide 
amounts in two independent samples. 

Models that define the variance as a function of the mean have been stud- 
ied intensively in the framework of heteroscedastic regression, where differ- 
ent estimation techniques have been suggested [e.g., Davidian and Carroll 
(1987)]. However, in MS, the variance function depends on the unknown 
peptide amount and the estimation problem is much more involved. A sim- 
ilar problem arises in certain immunoassay studies where few or no stan- 
dard concentrations are available [Raab (1981), Sadler and Smith (1986), 
O'Malley, Smith and Sadler (2008)], and in the evolving area of microarray 
mRNA expression analysis [Carroll and Wang (2008), Wang, Ma and Car- 
roll (2009), Fan, Feng and Niu (2010)]. Unlike immunoassay and microarray, 
labeled MS data can contain as few as two measurements of each peptide 
relative quantity and, therefore, analysis requires a special experiment for 
estimating the variance function. 

In a previous study, Zhang et al. (2010) suggest a novel controlled exper- 
iment for the estimation of the variance function in the iTRAQ protocol. 
Using the standard workflow, they held the sample constant by generating 
pairs of iTRAQ intensities from identical biological samples. Under this ex- 
periment, it is reasonable to assume that the iTRAQ labels are interchange- 
able in their error characteristics since labeled samples are mixed before 
processing and, therefore, the difference within pairs of such measurements 
are entirely attributable to the measurement error of the instrument itself. 
Two controlled experiments were conducted in January and March of 2009 
with the sole goal of estimating the variance function [Zhang et al. (2010)] . 
A separate experiment was conducted to explore differences observed be- 
tween wild-type and cancer cells expressing distinct mutations of the same 
oncogenic kinase (FLT3); this experiment relied on the variance function es- 
timated in the controlled study. The motivation for our current study arises 
from these past experiments, and below we describe in detail the mathemat- 
ical problem, suggest statistical methods to tackle it and apply them to the 
three data sets mentioned above. 

In Section 2 we present a naive method for estimating the variance func- 
tion employed by Zhang et al. and explore its validity. We show that the 
method works well when the error terms are typically small, as is the case 
in the instrument explored by Zhang et al. (2010), but may yield biased 
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estimators for the variance function in other situations. We then suggest an 
alternative mixture model approach for estimating the variance function and 
prove its consistency. In Section 3 we use the estimated variance function 
for interval estimation of the ratio of peptides across two different biological 
conditions, and we apply the method to iTRAQ-based MS analysis that is in- 
tended to decipher the oncogenic potential of two different clinically relevant 
mutations in the FLT3 receptor tyrosine kinase. Section 4 discusses testing 
of the hypothesis that the amounts of peptides in two biological samples are 
equivalent. The properties of the estimation approaches are investigated in 
Section 5 by simulation. Section 6 completes the paper with some remarks. 

2. Estimation of the variance function. 

2.1. The model. Consider a controlled iTRAQ experiment that quanti- 
fies N peptides. Let 1^1,1^2 denote two measures of intensity of peptide i 
(on a logarithmic scale) having an unknown mean //j. Assume that (l^i,li2) 
(i = 1, . . . , A^) are independent following the model: 

(1) Yij ~ N{ni,h{9, i^ii)), j = 1,2, independent, 

where is a vector of unknown variance parameters, and h(9,fM) is a known 
positive function, such as the power function Bi/j,^^ or the exponential func- 
tion e^"*"^^. In problems where /ii,...,/iAr are known or are modeled by 
a small number of auxiliary variables, standard techniques for estimating 
heteroscedastic regression models apply [e.g., Davidian and Carroll (1987)]. 
However, this is not the case in MS data where /ii,...,/iAr are unknown 
nuisance parameters. 

Klawonn, Hundertmark and Jansch (2006) and Hundertmark et al. (2009) 
developed an EM algorithm that maximizes over /xi , . . . , /xtv and 6 the likeli- 
hood corresponding to N independent observations, {yii,yi2), from model (1): 

(2) n{2vr/.(0,/..)} exp| ^^^^^-^ j. 

However, since the number of nuisance parameters increases with sample 
size, the maximum likelihood approach may provide biased estimators. The 
bias is readily seen in the classical example by Neyman and Scott (1948) 
of the one-parameter homoscedastic model h{9,fi) = 9. The maximum like- 
lihood estimator for 9 under this model is N~^ ^^^(^i ~ ^2)^/4 having ex- 
pectation 9/2, hence it converges to half the true variance. Thus, although 
the model is parametric, alternatives to the maximum likelihood technique 
should be employed. 

2.2. The MACL approach of Sadler and Smith. Motivated by immunoas- 
say data, Raab (1981) suggests to modify the likelihood (2) by multiplying 
the contribution of each of the pairs (1^1,1^2) by h^''^{9,fii), and then to 
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maximize the modified likelihood with respect to /ii,...,/X7v and 9. Raab 
shows by simulation that the standard maximum likelihood estimator is bi- 
ased, but the modified estimator performs reasonably well. Raab's method is 
computer intensive, as it requires the estimation of all nuisance parameters. 
Sadler and Smith (1986) estimate 9 by maximizing Raab's modified like- 
lihood at the point /ij = % := {Yn + Yi2)/2 {i = 1,. . . ,N): 

^ 1 

(3) 9 = argm^xll-—j==^exp{-Sf/2hi9,Y,)}, 

where Sf := {Ya - %? + (^,2 - %? = {Yn - Y^f jl. This approach, called 
maximum approximate conditional likelihood (MACL) by Sadler and Smith, 
reduces the estimation task to solving a small set of nonlinear equations. The 
resulting estimating equations under the normal model are therefore 

For example, for the model /i(0, ii) = exp{9i + 92fJ-), the estimating equations 
reduce to 

N 

(4) l-N-'Yl ^i exp(-ei - 92Y,) = 0, 

N N 

(5) N-'J2^i - N~'J2^iSfexp{-9i ~ 92?,) = 0, 

and can be easily solved by standard optimization algorithms using, for ex- 
ample, the R function optim [R Development Core Team (2011)]. As pointed 
out by Sadler and Smith, the solution for (3) can be obtained by an iterative 
reweighted least squares algorithm. 

We applied the MACL approach to the January (A^ = 2144) and March 
(A^ = 2174) experiments described in Section 1 using the functional form 
h{9,fi) = exp(6'i + 6*2/^) and obtained the estimates 9 = (4.89, —0.935) and 
9 = (4.89,-0.925), respectively.^ The similarity of the two estimated vari- 
ance functions is remarkable, suggesting that the between-study variability 
is small. This is a very important finding, as the variance function can be 
estimated in a control experiment and be applied to data obtained in inde- 
pendent experiments on the same instrument. We finally pooled the data 
together and obtained an overall MACL estimate of = (4.86, —0.927). 

Neither Raab nor Sadler and Smith provide sound theoretical justification 
for their methods, but explore them in several special relevant cases. In 



^Data and R codes can be accessed as project syn310406 on the Sage Bionetworks 
Synapse system (http://synapse.sagebase.org). 
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general, the expectation of the estimating equations differs from 0, hence, 
the estimators of the variance function are, in general, inconsistent. This is 
shown in Appendix A and is illustrated by simulation in Section 5. However, 
Appendix A suggests that the bias is small when the variance (for all /Xj) 
is small, because in such circumstances 1^ is a good estimator for /ij, even 
though it is based on only two observations. Recently, Wang, Ma and Carroll 
(2009) studied variance functions for microarray experiments and showed 
a similar inconsistency problem of estimators obtained by the method of 
moments. 

2.3. A mixture model. A possible strategy to deal with the inconsistency 
of the MACL approach is to impose additional reasonable assumptions on 
the nuisance parameters. We consider the model 

Yijlfiir^ N{fii,h{9,fii)), J = 1,2, independent, 

(6) 

fii r^ Go, i = 1,. . . ,N, independent, 

where (i) the support of Go is in the segment [a, b], that is, P{a < /Uj < 6) = 1, 
(ii) the variance is bounded, that is, a < h{9, n) < /3 for all // in the support of 
Go for some < q < /3 < oo, (iii) h{6,fi) is continuous on [a,b] and identifies 
6, that is, knowing h{9,n) on the support of Gq implies knowledge of 6. 

The assumptions on h are satisfied by most practical models. The rea- 
son for bounding Go and the choice of the values a and b are discussed in 
Section 3. To see the importance of the identifiability assumption, consider 
the model h{0,^) = exp(6'i + ^2/^) and a degenerate Go that assigns all the 
mass to some /Uq. In such a model, exp(0i + ^2/^0) is identifiable, but the 
pair (^1,^2) is not. 

Theorem 1. Under model (6) and the assumptions following it, the 
maximum likelihood estimator of {6, Gq) is consistent. 

The proof, which is sketched in Appendix B, is based on the seminal paper 
by Kiefer and Wolfowitz (1956) who prove the consistency of the maximum 
likelihood estimator in mixture models such as (6). A recent application of 
mixture models for variance estimation in microarray analysis can be found 
in Wang, Ma and Carroll (2009), though they suggest a different estimation 
strategy for 9. 

Several algorithms for deriving the maximum likelihood estimator of a 
mixture model have been suggested in the literature [see, e.g., Bohning 
(1999)]. Here we estimate 9 and Gq by employing the EM algorithm, treat- 
ing {Yii,Yi2,fJ.i) as the complete data on the ith peptide. One strategy for 
estimating Gq is to restrict the search to distributions supported on a fine 
grid and to find the maximum likelihood among these distributions. In the 
current problem, the variance becomes small for large values of fx and using 
a simple grid may lead to data points which are too far (in terms of stan- 
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dard deviations) from all support points. We found that defining the support 
points of Go as a function of the variance performed better than using a sim- 
ple grid. Thus, we first obtained an initial estimate of the variance function 
using, for example, the MACL approach, and then restricted the distance 
between two support points to be at most d standard deviations according 
to the estimated function. Specifically, let 6 be an initial estimate for 9, 
and define the maximal support point to be fij = b, the second to maximal 



to be hj-i = fij — d\/ h{6,^j) and recursively fJ-j-i = jJij — dy h{9,iJ,j). The 
selected points of support depend on the initial estimate of the variance 
function and can be updated as part of the algorithm, though in our ex- 
perience, this update made no significant improvement when using d = 1/4. 
Once the support points for Go are determined, the EM algorithm is applied 
to estimate 6 and Go. The algorithm is quite standard and its description 
is detailed in Appendix C. 

We fit the following three forms for h{9,fi) that have been suggested 
in the literature^: exp{9i) fj.^'^ , exp(^i -|- 02M) ^iid exp(0i -|- ^2/^) + exp(6'3). 
The right panel of Figure 1 presents estimates of these three functions on 
a scatter diagram of Yi = {Yn + Yi2)/2 against Sf = (Ya - Yi2f/2. The X- 
axis is a naive estimate of the mean of each pair, /Xj, and the Y-axis is a 
naive estimate of the variance. All three models give similar results in most 
of the range with some deviation for very small values of /i. A loess fit, 
presented in the figure by the dotted line, is quite close to the functional 
form h{9,fi) = exp{9i + ^2A*) originally suggested by Zhang et al. (2010). 
This latter model is used throughout this paper. 

Applying the EM approach using the functional form h{9,iJ,) = exp{9i + 
92fi), we obtained the estimates 9 = (4.91,-0.929) and 9 = (4.96,-0.944) 
based on the January and March experiments, respectively. The estimates 
are very similar to the estimates obtained by the MACL approach. We then 
pooled the data from the two experiments together and obtained our final 
estimate, 9 = (4.84, —0.927). The corresponding estimate of Go is displayed 
in the supplemental article [Mandel et al. (2013)]. Initial values for the EM 
algorithm were obtained by the MACL approach. Starting the algorithm 
from different points resulted in essentially the same estimate (details are 
provided in the supplementary materials). 

3. Confidence intervals. Having estimated the variance function, confi- 
dence intervals for various parameters of interest can be constructed based 
on data from a new experiment that compares different biological samples. 
We construct frequentist confidence intervals that are based on the estimate 



■^Data and R codes can be accessed as project syn310406 on the Sage Bionetworks 
Synapse system (http://synapse.sagebase.org). 
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6 oi 6 under model (6), but do not use the estimate of the mixing distribution 
Gq. We take this approach since the variance function is a stable character- 
istic of the mass spectrometer, whereas Go depends on the biological sample 
and may differ from sample to sample. 

3.1. Confidence intervals for //. It is of interest to attach a measure of 
uncertainty to the observed intensity or to report a range rather than one 
value for each peptide. This section discusses the construction of 1 — a confi- 
dence intervals for the peptide amount, /i, based on one observation Y from 
the model K ~ N{^,h{6,^)); the next section deals with the construction 
of confidence sets for the relative abundance of a peptide in two different 
samples. 

Similar to the MACL approach for estimation, construction of confidence 
intervals can be simplified by plugging an estimate of fi in h{6,fi). Thus, 



a naive confidence interval is constructed by y it zi„a/2y /i(^,5^), where 
Za is the a quantile of the standard normal distribution. This method is 
expected to perform reasonably well only in cases were the variance is small 
and changes slowly with //. 

An exact 1 — a confidence set for /x can be constructed using the piv- 
otal quantity grip-) = (Y — fi)'^ / h{9 , n) . This interval is defined as Ca = 
{fi'-gyifJ-) < Xi,i-q}) where Xdf,a is the a quantile of the chi-squared distri- 
bution with df degrees of freedom. 

For the model h{6,fi) = exp{6i + 02fJ.), this set can be easily found by a 
bisection search using the following observations: 

(1) gviY) = is a local minimum. 

(2) For ^2 < 0, /i* = Y + 2/62 is a local maximum of 5fy(/i), with gyifJ-*) = 
402-2e-(2+ei+^2y). 

We thus obtain the following properties of the confidence procedure: 

• If gvifJ-*) ^ Xi-a, then the confidence set is a one-sided interval of the 
form (— oo,ri). 

• If gvilJ*) > Xi-cn then the confidence set is a union of two intervals, 
(-00, n) U (^2,^2), where ri< fi* <l2<Y < r2. 

This nonstandard shape of the confidence set refiects the fact that a realiza- 
tion y is likely either when fj, is close to y or when fx is much smaller than 
y and the variance of the measurement is very large. Often, the range of 
^ is a priori bounded so that values smaller than ri do not belong to the 
parameter space, and the confidence set is always an interval. 

Since the parameter 6 is unknown, a consistent estimator based on the 
mixture approach is plugged in to generate confidence intervals with an 
approximate level 1 — a. 
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3.2. Confidence intervals for fj,i — fi2- Let Yi ~ iV(/ii,/i(0,/ii)) and I2 ~ 
N{fi2,h{6,H2)) be the log intensities of the same peptide obtained by the 
iTRAQ protocol under two different conditions, and assume that Yi and I2 
are independent conditionally on ni and fj,2- We consider the construction 
of a confidence set for fii — /U2, which is the parameter of primary interest 
in many quantitative MS studies. 

As in the one-parameter case, a naive confidence interval for fii — ^2 
can be obtained by plugging Yi and Y2 into the variance term: Yi — I2 i 



^i~a/2\J ^{^1 Yi) + h{9, Y2). However, such intervals may be anti-conservative, 
as the estimator Yi for /ij is inconsistent, and hence so is the estimator 

h{e,Yi)iovh{e,^ii). 

A direct and a relatively simple way of calculating conservative intervals 
is by Bonferroni correction, that is, by first constructing 1 — a/2 intervals for 
/ii and for /i2, as described in the previous section, and then calculating the 
minimum and maximum differences of the two intervals. However, a more 
direct construction uses the reparametrization z^i = ^1 — ^2 and V2 = ^1 + ^2- 

Let 

Yi - Y2 - VI 



and 



C/yj-Yj (1/1,1^2) 



9Yi+Y2{vi,^2) 



y/h{9, (y2 - yi)/2) + h{e, {V2 + i^i)/2) 

Yi + Y2 - V2 

^h(B, {V2 - T^l)m + h{e, {V2 + V{)I1) ' 



then {gY-^-Y^{y\^V2)^gY^+Y2(y\^V2)) has a bivariate standard normal distri- 
bution with correlation 

p{yi.V2) = {hie, {U2 + yi)/2) - hie, {V2 - z/i)/2)} 

/{h{e, {v2 - v^)i2) + h(e, (v2 + z^i)/2)}. 

Thus, 

9Yiy2{v\,V2) = {gYi-~Y2{vU^2),gYi+Y2{vi,^2)) 

1 P{t^I,V2)\ f gYi~Y2{vU^2) 

.P(l^l,l^2) 1 J V 5^1+^2(1^1, 1^2), 

has a Xj2) distribution and can serve as a pivot for constructing confidence 
regions. Specifically, 

(7) Cy^,Y2{^Ii^2) = {{^1,^2) ■■ gYi,Y2{^l^^2) < X2,l~a} 

is an exact 1 — a confidence region for {1^1,1^2) and, therefore, {vi : (j^i, 2^2) £ 
C'Yi,Y2i'^iy'^2), —00 <V2 < oo} is a Conservative confidence set for vi. 

3.3. An application to the iTRAQ protocol. In this section we analyze 
data from the two control experiments mentioned in Section 1. As the two 
biological samples in these experiments were identical, nn — pi2 = for all 
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Fig. 2. 95% confidence regions for the pair (v\,V2) constructed for two data points: (8, 9) 
(left) and (7.5, 8) (right). 



peptides i = l, . . . ,N.ln order to evaluate the performance of our confidence 
intervals, we used the parameters' estimates from one experiment to con- 
struct 95% level confidence intervals for the difference of peptide abundances 
in the other experiment, and calculated the proportion of intervals that did 
not cover 0, the true difference. 

As in the uni-parameter case, the confidence set for (1^1,1^2) is not neces- 
sarily a connected set and the resulting confidence set for ui is not always 
an interval. Figure 2 demonstrates the shape of the confidence region for 
the pairs (^1,12) = (8,9) and (^1,^2) = (7.5,8). The X and Y axes are, re- 
spectively, z^i and z^2) and the confidence sets are the shaded areas. The 
figure reveals that without restricting the parameter space of /i (hence the 
parameter space of z^i), the confidence interval for z^i comprises all of the 
real line and is noninformative. We therefore assumed that fi £ [7.3,13.9], 
where the limits were determined from typically observed values as well as 
the expected range of the intensity values. This decision restricts the values 
of i^i and 1^2 to the large parallelograms depicted in Figure 2, and enables 
the construction of informative confidence intervals. The smaller parallelo- 
grams represent confidence sets obtained by Bonferroni correction applied 
to univariate confidence intervals for /xi and for fi2- 

Using estimates from the January data, we calculated confidence intervals 
for vi for peptides in the March experiment by (7), Bonferroni correction 
and the naive approach, and found that, respectively, 18 (0.8%), 1 (0.05%) 
and 86 (3.96%) of the 2174 intervals did not include the true parameter 
ui = 0. The corresponding numbers for the 2144 peptides in the January 
experiments using estimates from the March data are 47 (2.19%), 8 (0.4%) 
and 106 (4.94%). This exercise suggests that interval (7) is better than the 
Bonferroni interval, but is still conservative. The naive approach performs 



VARIANCE FUNCTION ESTIMATION IN QUANTITATIVE MS 11 

Table 1 

95% confidence intervals for the ratio of phosphopeptide quantity across two experimental 

conditions using the conservative and the naive approaches 



Phosphopeptide 


FLT3-D835Y 


FLT3-ITD 


CI 


CI naive 


VLPQDKEpYYK 


10.21 


10.78 


(0.41,0.76) 


(0.44,0.72) 


GQESEpYGNITYPPAVR 


13.62 


11.89 


(5.05,6.36) 


(5.11,6.22) 


HKEEVpYENVHSK 


11.19 


9.92 


(2.66,5.05) 


(2.76,4.59) 


pYKNILPFDHSR 


10.83 


9.80 


(2.03,4.10) 


(2.12,3.69) 


AVDGpYVKPQIK 


11.45 


13.36 


(0.13,0.17) 


(0.13,0.17) 



surprisingly well in our experiment, but, in general, its theoretical coverage 
probability is not controlled. Further research is needed to understand this 
phenomenon. 

3.4. An application to cancer phosphoproteomics. Zhang et al. (2010) 
used iTRAQ labeling to study a key modification to proteins called phos- 
phorylation, which is important in cell signaling, and is often deregulated in 
cancer cells. Their study focused on aberrant signaling arising from onco- 
genic FLT3 mutations in acute myeloid leukemia. In particular, they moni- 
tored a key, subcomponent of signal transduction, namely, protein tyrosine 
phosphorylation, in order to obtain a global understanding of the onco- 
genic potential of two clinically identified FLT3 mutants (FLT3-D835Y and 
FLT3-ITD). Both FLT3 mutants induce constitutive signaling even in the 
absence of proper external cues, which results in uncontrolled cell prolifera- 
tion, a hallmark of cancer development [Blume-Jensen and Hunter (2001)]. 
Since receptor tyrosine kinases are often constitutively active in cancer cells, 
it is not surprising that a large proportion of downstream phosphorylation 
events diverge from a 1:1 value (measured relative to a control cell line), 
effectively eliminating the possibility of deriving a variance function from 
the experimental data themselves. 

Based on the mixture model results of the pooled control experiments, we 
calculated 95% confidence intervals for the ratio of phosphopeptides men- 
tioned explicitly in the figures of Zhang et al.; these are reported in Ta- 
ble 1. For comparison, we calculated confidence intervals based on a naive 
approach mentioned in Section 3.2. Using the intensity-dependent variance 
function described in this paper, as opposed to a constant cutoff point often 
used in the literature, subtle changes in phosphorylation levels can be found 
(e.g., peptide VLPQDKEpYYK). Such ratios can be considered statistically 
significant despite being smaller than the typical cutoff of 2:1. This in turn 
suggests that experiments can be designed to explore phosphorylation under 
physiological conditions, unlike many current experimental designs where ar- 
tificial or extreme environments are used in order to amplify the observed 
ratios. 
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4. Hypothesis testing. 

4.1. Calculation ofp-values. For a given peptide, consider testing the hy- 
pothesis Hq : Hi = fi2 = fJ- versus the two-sided alternative i/i : /ii 7^ /i2- The- 
oreticahy, this can be done by inverting the confidence intervals discussed 
in the previous section. However, for calculating p-values, this inversion is 
computationally difficult and the current section explores an alternative ap- 
proach. 

As before, let Yi ~ N{fii,h{6,fii)) and Y2 ~ ]^{l^2,h{6,fi2)) be indepen- 
dent, then Yi — I2 ~ -^(^1 — lJ'2,h{6,fii) + h{6,H2)) and a reasonable test 
may compare 

(8) (^1-^^)' 

to the chi-squared distribution with one degree of freedom. However, (8) 
contains the unknown parameters 9 and /u, and hence is not a legitimate test 
statistic. As in the construction of confidence intervals, a naive p- value can 
be calculated by replacing the denominator of (8) with 2h{6,{Yi +Y2)/2). 
Although it is reasonable to replace 9 with its consistent estimator 6, the 
average (li + l2)/2 is an inconsistent estimator for fi, possibly leading to 
an anti-conservative p-value. 

As in other statistical problems involving nuisance parameters, an asymp- 
totically conservative p-value is obtained by 

sup P{Z^>{yi-y2f/2h{9,fi)), 

where yi and 7/2 are the realized values and Z"^ is a xfi) random variable. 
This approach is simple and easy to implement but results in an overly 
conservative p-value. We therefore suggest to employ the approach of Berger 
and Boos (1994), where the p-value is calculated by maximization over a 
confidence interval for the nuisance parameter. Specifically, let Cphe al — /3 
level confidence interval for //, obtained in a way similar to that presented 
in Section 3.1, then we define our p-value as 

(9) snpP{Z^>{yi-y2f/2h{9,fi)) + (3. 

The choice of /3 depends on the context. If a univariate hypothesis is tested 
with a significant level of 5%, then /3 = 0.001 is usually a good choice. How- 
ever, if a Bonferroni correction is needed, then /3 must be much smaller, as 
p-value > /5 by construction. 

4.2. Application to the iTRAQ data, p-values for testing no difference 
of peptide amounts were calculated for all pairs in the pooled iTRAQ con- 
trol experiments described in Section 3.3. We used the variance function 
h{9,fi) = exp(6'i -|- 6*2/-^) with 9 estimated by the EM algorithm applied to 
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Fig. 3. Scatterplots and histograms of p-values for the pooled control experiments ob- 
tained by different approaches. Scatterplots are in log scale with dotted lines indicating 
significant level of 5% with and without Bonferroni correction. 



the pooled data. Three approaches were compared: the naive approach that 
replaces // in h{6,^) with Y = (Yi + 5^2)/2, a conservative simple approach 
that replaces /z with a, and the approach of Berger and Boos based on (9) 
with /? = 10-6. 

Figure 3 presents scatterplots and histograms of the p-values. Scatter- 
plots are depicted in logarithmic scale with dotted lines indicating 0.05 and 
O.OS/A'^ significant levels. Recall that the null hypothesis holds in these exper- 
iments. Two peptides had naive p-values smaller than the Bonferroni cutoff 
value, for one of them the Berger and Boos p-value was also under the cutoff 
level. For about 40% of the experiments the Berger and Boos p-value was 
similar to the conservative p-value, but for the other peptides, the Berger 
and Boos method gives considerable smaller values. These are valid p-values 
and certainly worth the additional computation effort. The distribution of 
p-values under the naive approach is very close to the uniform distribution, 
whereas the distributions under the other two approaches are stochastically 
larger. This and the simulation presented in the next section suggest that 
the naive approach for testing differences of peptide amounts may be only 
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slightly anti-conservative, though a more extensive study is required before 
the naive approach can be recommended. 

As a second application, we calculated p-values for the data described 
in Section 3.4 that contrast two mutants found often in certain types of 
cancer. There are A^ = 205 peptides in the data, the p- values of 141, 123 and 
9 of them were smaller then 0.05/A according to the naive, the Berger and 
Boos and the conservative approach, respectively. These peptides are then 
prioritized for in-depth functional characterization. 

5. Simulation. 

5.1. Estimating the variance function. The performances of the two es- 
timation approaches, the MACL and the mixture model, were tested by 
simulation under various conditions. The first set of simulations aimed at 
testing the performance of the MACL approach. For each of the scenarios 
described below, we generated values for the nuisance parameters and then 
generated independent pairs of observations from the corresponding normal 
distributions. We repeated this process 1000 times for different sample sizes 
{N = 200, 500, 1000 and 2000) and for two different variance functions of the 
form exp(^i + ^2/^): (^ii ^2) = (5, —1), which is similar to the values obtained 
in our data, and (^1,^2) = (5, —0.5), which reflects observations with a much 
larger variance. The following scenarios were considered: 

• Observations fixed: a set of /ij values was sampled from the observed l^'s 
(with replacement) and the same values were used in all 1000 replications. 

• Observations random: a different set of /Uj values was sampled from the 
observed 1^'s (with replacement) for each simulation. 

• U(8,12): the fii values were generated from the continuous uniform dis- 
tribution over (8,12). 

• U{8,9, . . . , 12}: the fii values were generated from the discrete uniform 
distribution over {8, 9, . . . , 12}. 

The results of the simulation are summarized in Table 2 and in more 
details in the supplemental article [Mandel et al. (2013)]. There seems to 
be almost no difference between the scenarios considered (see Table 1 of the 
supplementary materials), and this suggests that the approach is insensitive 
to modest changes in the distribution of the nuisance parameters. Both the 
variance and the bias decrease with sample size for the model (^1,^2) = 
(5,-1), and the overall performance of the MACL approach for this case is 
satisfactory. However, the MACL estimators are biased for the case (0i, ^2) = 
(5,-0.5), and the bias did not decrease with sample size (Table 2). Thus, 
unless the variance is very small, the approach is problematic and is not 
recommended. 
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Table 2 
Simulation results for the MACL method under the observations fixed scenario 







01 






02 




N 


6*1 


bias 


std 


6>2 


bias 


std 


200 


5 


-0.061 


0.830 


-1 


0.006 


0.081 


500 


5 


-0.051 


0.508 


-1 


0.005 


0.049 


1000 


5 


-0.020 


0.360 


-1 


0.002 


0.036 


2000 


5 


-0.002 


0.251 


-1 


0.000 


0.025 


200 


5 


-1.238 


0.759 


-0.5 


0.127 


0.074 


500 


5 


-1.175 


0.484 


-0.5 


0.121 


0.047 


1000 


5 


-1.156 


0.319 


-0.5 


0.119 


0.031 


2000 


5 


-1.164 


0.238 


-0.5 


0.120 


0.023 



In order to test the mixture model approach, we generated /ij by the ob- 
servations fixed scenario under the two variance functions described above. 
For each sample size, we simulated 200 data sets and calculated the em- 
pirical biases and standard deviations. The results are listed in Table 3. As 
expected, the bias and variance of the estimators decrease with sample size 
for both models. 

Figures 4 and 5 display the performance of the estimators for the variance, 
that is, the performance of exp(^i + 62IJ) as a function of ^. The gray lines 
are estimated variance functions from 200 simulated data sets and the 
true variance function is depicted in black. The figures demonstrate that 
the mixture model approach is as good as the MACL in the low variance 
case and performs better in the large variance scenario, especially when the 
sample size is large. 

5.2. Confidence intervals for /x. Intervals for /i based on one observation 
use the pivot {Y — n)'^/h{6, fj) and are exact. Here we study the performance 

Table 3 
Simulation results for the mixture model method under the observations fixed scenario 







01 






02 




N 


ei 


bias 


std 


02 


bias 


std 


200 


5 


0.580 


0.914 


-1 


-0.071 


0.088 


500 


5 


0.423 


0.534 


-1 


-0.049 


0.052 


1000 


5 


0.274 


0.328 


-1 


-0.030 


0.032 


2000 


5 


0.173 


0.227 


-1 


-0.019 


0.022 


200 


5 


0.070 


1.026 


-0.5 


-0.009 


0.099 


500 


5 


0.019 


0.608 


-0.5 


-0.003 


0.059 


1000 


5 


0.027 


0.397 


-0.5 


-0.003 


0.039 


2000 


5 


-0.013 


0.291 


-0.5 


0.001 


0.028 
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Fig. 4. Estimated variance functions (gray) and the true variance function (black) ob- 
tained in 200 simulated data sets for the case {61,82) — (5,-1). The figures on the left 
show the results of the MACL approach and those on the right are for the mixture model 
approach. The simulated data sample sizes are, from top to bottom, 200, 500, 1000 and 
2000. 
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Fig. 5. Estimated variance functions (gray) and the true variance function (black) ob- 
tained in 200 simulated data sets for the case {61,82) = (5, —0.5). The figures on the left 
show the results of the MACL approach and those on the right are for the mixture model 
approach. The simulated data sample sizes are, from top to bottom, 200, 500, 1000 and 
2000. 
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Fig. 6. Coverage probability of the naive confidence interval for /j, for confidence levels 
0.99, 0.95 and 0.90. Solid lines: h{6, n) = exp{5~ fi), dashed lines: h{6,ij,) =exp{5 — 0.5fi). 



of the corresponding 1 — a naive confidence intervals that replace fi with Y 
in the variance function h{6,fj,); see Section 3.1. 

Figure 6 presents the coverage probability of the naive intervals as a 
function of the mean (// = 7, 7.1, . . . , 14), the coverage probability (1 — a = 
0.9, 0.95, 0.99) and the variance function [exp(5 — fi), exp(5 — 0.5/x)]. For each 
fi, a and variance function h{6,fi), we estimated the coverage probability by 
simulating 100,000 replications from the model N{fi,h{9,fi)), constructing 
naive confidence intervals of level 1 — a, and calculating the proportion of 
intervals covering //. The performance depends on the variance at //, where 
the true coverage for small /x (i.e., a large variance) could be much lower 
than the aimed coverage, especially for large values of 1 — a and for the 
model h{9,fi) = exp(5 — 0.5/i) represented by dashed lines. 

5.3. p-values. A third simulation study was conducted with the aim 
of understanding the properties of p-values obtained by the different ap- 
proaches. For a selected set of values for fi, we generated 10,000 inde- 
pendent pairs of observations (Yi,Y2) such that Yi ~ N{n,e^~^ '^^) and 
12 ~ ^(/Wfc,e^^'*'^^^''), independently, where Hk = ii + k x \/e^i+^, that is, 
Yi and Y2 are centered about k standard deviations apart. The parameter 
k ranged from to 3, and for (^1,^2) we studied the values (5,-1) and 
(5, —0.5). We used the value j3 = 10~^ for the Berger and Boos p-values. 

Figures 4 and 5 of the supplemental article [Mandel et al. (2013)] present 
the proportions of p-values that were smaller than 0.05 as a function of 
11, k and 62- For the case 82 = —1 (supplemental article. Figure 5), the 
naive approach is only slightly anti-conservative and only for small values 
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of fi, and its power is much larger than that of the other methods. The 
Berger and Boos approach works reasonably well only for very small and 
very large values of /x; the conservative approach is useless. For the case 62 = 
—0.5 (supplemental article, Figure 4), the naive approach is anti-conservative 
for small values of // with a significant level of up to 0.08 instead of the 
declared level of 0.05. However, it is the only method that has a useful 
power function. 

6. Discussion. This paper presents a protocol for estimating the variance 
function of a mass spectrometer when used for relative quantification in pro- 
teomic applications. Using two sets of data collected three months apart, we 
found that the variance function is stable over time. However, we expect to 
find different variance functions in different instruments and, hence, each lab 
should estimate the variance parameters of each spectrometer independently 
using the protocol described here, and update them periodically. More im- 
portantly, the variance estimated here corresponds only to the variance of 
the instrument itself, and does not include error terms corresponding to 
the natural variability of biological samples or the processing required to 
solubilize proteins from cells and tissues, generate peptides, etc. 

Two inference approaches are considered. The first estimates the nuisance 
parameters by simple averages and then maximizes a target function; the 
second approach assumes a mixture model and estimates the variance func- 
tion by maximum likelihood. When the variance is small and changes slowly 
as a function of the mean, as is the case in the data we analyzed, an average 
of two iTRAQ reporter ions per peptide provides a reasonable estimate for 
the unknown /j. and the first approach gives good results. However, it may 
yield highly biased estimators in other scenarios, as demonstrated by simu- 
lations, and we therefore recommend the routine use of the mixture model 
approach that has a sound theoretical justification. 

When using the estimated variance function for statistical inference on 
one parameter /x, exact methods for constructing confidence intervals are 
available and are much more appropriate than intervals constructed by the 
anti-conservative naive approach. On the other hand, for inference on the 
difference or the ratio of peptide abundance measured across two biological 
conditions, the naive approach performs quite well, yielding a significant 
level only slightly larger than the aimed one. 

The iTRAQ protocol is somewhat more complicated than presented here, 
as it involves preprocessing of the spectral data to correct for differences 
in total protein amount, iTRAQ label purity and instrument-specific pa- 
rameters. Moreover, iTRAQ is known to suffer from contamination due to 
co-eluting chromatographic peaks that also share similar precursor masses 
(i.e., peptides which have a similar m/z value and retention time). Theoret- 
ically, these factors may induce dependence between measurements that is 
ignored in the current analysis. 
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In order to test the underlying normal assumption, Zhang et al. produced 
a q-q plot of {Yn — Yi2)/y h{9,Yi) that showed a very good fit. Although 

the graph is suggestive, it relies on Yi as an estimate for the nuisance pa- 
rameter fii, hence it does not have a theoretical support. A formal approach 
we intend to explore requires at least three observations for each peptide. 
A simple transformation of each of the triplets results in variables that, un- 
der the normal model, have a Cauchy distribution, and a q-q plot or formal 
goodness-of-fit tests can be easily employed. This approach requires iTRAQ 
data from a control experiment that yields more than two independent and 
identically distributed measures for each peptide. We intend to conduct such 
an experiment in the future. 

The conservative results of the exercise conducted in Section 3.2 are par- 
tially due to the inherently conservative construction of the intervals, but 
may also be a result of the special parametric shape for the variance func- 
tion we considered. Our chosen parametric model is very simple, enabling 
the implementation of simple algorithms, and is supported quite well by the 
data. A nonparametric method has been recently suggested for the analysis 
of microarray data [Carroll and Wang (2008)]; the possibility of adopting 
this method to MS data and of using it for goodness-of-fit testing should be 
further explored. 

APPENDIX A: BIAS OF THE ESTIMATING EQUATIONS 
Recall that Yi ~ N{ni, le^i+^2M») and Sf ~ e'^^+^^^''xfl), and that Yi and 

Sf are independent. Straightforward calculations show that E{Sf) = e^^ '^^^ 
and ^{exp(-6'i - 62?,)} = exp(-6li - e2^ii + le^e^^+^^^'^) and, therefore, 

i?|iV-i|;52exp(-0i-e2F.)|=iV-i|;expQ02V^+^^^») 

so that the expectation of the first equation (4) differs from zero, unless 
62 = 0, that is, the homogeneous model of Neyman and Scott (1948). 

For the second equation, we have EYi exp(— ^2^^) = — -^^E exp{—tYi)\t=g2 = 
_^exp{-t/i, + t2ie^i+^-/^'}|t=e2 = ifii - i02e^i+^2'^Oexp{-02Aii + Ofj x 
eei+e2/.,|^ g^ ^^^^ E{YiSf expi-Oi - 62%)} = i^l^ - ^^2e^^+^'^Oexp(i0i x 
e^i+^2^*), and the expectation of the left-hand side of (5) is 

AT-i f;/., - N''Y.U - W'^^'^^^) exp(i^ie^^+^^^») , 

which again differs from for ^2 7^ 0. 

In general, the bias will be small if e^"*"^^' is small for all i, which means 
that the variance of the measurements is small and, hence, the local averages 
are good estimators for the unknown ^j parameters. 
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APPENDIX B: CONSISTENCY OF THE MLE 
A generic sample point isy = (yi, 2/2)) where, by conditional independence, 

(10) /(,;%) = {2vr/.(^,^)}-exp{-i^l^i|±i|^i^ 

The marginal density ofY={Yi,Y2) is g{y;e,Go) = fj {y; 6\t) dGo{t) . 

The proof of consistency is based on the result of Kiefer and Wolfowitz 
(1956) (KW hereafter); the metric we use below is given in KW equa- 
tion (2.2). 

We complete the parameter space of {6, Gq) by including all proper dis- 
tributions functions with support in [a,b]. 

Assumption 1 of KW trivially holds with respect to the Lebesgue measure. 
Next, note that f{y;9\fi), and hence g{y;6,Go), is bounded above by a~^ . 
Let {9i,Gi) — >■ {9*,G*), where {9*,G*) is in the complete parameter space. 
In order to verify Assumption 2 of KW, we need to show that g^y; 9i,Gi) — )• 
g{y;9*,G*). We have 

\giy;9,,G,)-g{y;9*,G*)\ 



< 



f{y;9i\t)dG,{t)- J f{y;9*\t)dG*{t) 
{fiy;9,\t)-fiy;9*\t)}dG*it) 



+ 



f{y;9i\t)d{G,{t)-G*{t)) 



<J\f{y;9i\t)-fiy;9*\t)\dG*{t) + a-^Jd\G,{t)-G*{t)\. 

The first term vanishes by the Dominated Convergence theorem and the 
second vanishes by the convergence of Gi to G* . 

For verifying Assumption 3 of KW, define 'm{y; 9*,G* , p) = sup g{y; 9, G), 
where the supremum is taken over all {9,G) such that \9 — 9*\ + \G — G*\ < p. 
We need to show that m is a measurable function of y for any p > and 
any {9*,G*) in the complete parameter space. This is true for the same 
arguments given by KW in their first example: g is for each y continuous in 
{9,G) and the parameter space is separable. To show this formally, define 
A{9*,G*,p,c) = {y:m{y;9*,G*,p) > c}, and let {{9i,Gi)} and {yj} be dense 
subsets in the parameter and sample space, respectively. Let B{y,r) and 
B{9,G,r) be balls of radius r around the corresponding points, then 

00 
A{9*,G*,p,c)=f]\jB{y„l/n), 

n=l j 

where the union is over {j : 3{9i,Gi) G B{9*,G* ,p) such that g{yj;9i,Gi) > c}. 
Assumption 4 of identification follows from Bruni and Koch [(1985), The- 
orem 1], that proves that Gq and h{9, p) are identifiable on the support of p. 
Assumption (iii) below equation (6) ensures identifiability of 9. 
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To verify Assumption 5, note that our assumptions on h{6,fi) guarantee 
that g{y; 6, G) is bounded above and below so that Elog{g{Y; 9, G)} > — oo, 
where the expectation is taken with respect to g{y;9Q,Go), the true density 

of y. 

APPENDIX C: AN EM ALGORITHM 

Let f{y;9\fi) be the bivariate normal density defined in (10), and let 
a < Hi < • • • , fj,j < b he fixed scalars (support points of Go). We approximate 
the likelihood of one observed pair by the following discrete mixture model: 

J 

(11) g{y;6,7v) = J2^jf{y-d\f,j), 

i=i 
where tv = (vri, . . . , vrj), ttj > and tti + • • • + ttj = 1. The unknown param- 
eters are the ttj's and 9. 

To construct the EM algorithm, consider a Multinomial variable A over 
1, . . . , J with a probability vector vr, and define {5i,. . . , 5j) by 6j = /{A = j}, 
where / is the indicator function. Let Y = (yi, . . . ,yj\f) be data on N pairs, 
then the complete log likelihood can be written as 

N J N J 

(12) £(7r,e;y) = j;j;5i,log{/(y,;e|^,)} + ^^5,,log(7r,), 

i=l j = l i=l j = l 

where 6ij is an indicator for the (unobserved) event {pair i has mean fij}. 

Denote by old the current estimates of the unknown parameters, then, 
using the Bayes formula, the expectation step reduces to estimating 

E-'''{6,,\Y) = E°'''{6ij\y,) = P"'\6ij = l\y,) 
7rf/(y,;0°i<i|/i,) ^^, 

— — Wi 



Note that X]7=i ^ij = 1 by definition, so the above formula can be interpreted 
as the current estimate of the probability that yi was generated by the 
distribution having mean fj,j. 

The maximization step is obtained by replacing 6ij in (12) with w°j and 
solving 

N J N J 

(13) maxJ]J]^°|'ilog{/(y.;0|^.,)} + J]J]u;°)'^log(7r,), 

which is done separately for vr and 9. For 9, the problem is of a nonpar ametric 
regression type and can be solved by reweighted least squares, similar to the 
MACL approach. The mixing probabilities are simply updated by 

1 ^ 
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SUPPLEMENTARY MATERIAL 

Web-based supplementary materials variance function estimation in quan- 
titative mass spectrometry with application to iTRAQ labeling 

(DOl: 10.1214/12- AOAS572SUPP; .pdf). Section A: Workflow of the iTRAQ 
technique. Section B: Estimate of Go- Section C: Sensitivity of the EM al- 
gorithm to initial values. Section D: Detailed simulation results. 
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